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ABSTRACT 

Cold fronts (CFs) are found in most galaxy clusters, as well as in some galaxies and 
groups of galaxies. We propose that some CFs are relics of merging between two 
shocks propagating in the same direction. Such shock mergers typically result in a 
quasi-spherical, factor sa 1.4 — 2.7 discontinuity in density and in temperature. These 
CFs may be found as far out as the virial shock, unlike what is expected in other CF 
formation models. 

As a demonstration of this effect, we use one dimensional simulations of clusters 
and show that shock induced cold fronts form when perturbations such as explosions 
or mergers occur near the cluster's centre. Perturbations at a cluster's core induce peri- 
odic merging between the virial shock and outgoing secondary shocks. These collisions 
yield a distinctive, concentric, geometric sequence of CFs which trace the expansion 
of the virial shock. 

Key words: galaxies: clusters: general - galaxies: haloes - X-rays: galaxies: clusters 
- shock waves 



1 INTRODUCTION 

Recent X-ray observations of galaxy clusters reveal various 
phenomena in the gaseous haloes of clusters, such as merg- 
ers, cavities, shocks and cold fronts (CFs). CFs are thought 
to be contact discontinuities, where the density and tem- 
perature jump, while the pressure remains continuous (up 
to pr ojection/resolution effects). Th ey are common in clus- 
ters (|Markevitch fc Vikhlininl l2007i ) and vary in morphol- 
ogy (they are often arcs but some linear or filamentary 
CFs exist), contrast (the density jump is scattered around 
a value of ~ 2), quiescence (some are in highly disturbed 
merging regions, and some in smooth, quiet locations), and 
orientation (some are arcs around the cluster centres, and 
some are radial, or spiral in nature). CFs have been postu- 
late d to originate from cold material stripped during merg- 
ers (jMarkevitch et all |2000| . for example) or slosh i ng of 
the intergalactic medium (IG M; iMarkevitch et all l200ll : 
lAscasibar fc MarkevitcfJ l200ot ). In some cases, a metallic- 
ity gradient is observed across the CF, indicating that it 
results from stripped gas, o r radial motions of gas (see 
IMarkevitch fc Vikhlininll2007l . and references therein). Shear 
is often found along CFs in relaxed cluster cores, implying 
nearly s onic bulk flow bene ath (towards the cluster's core) 
the CF (|Keshet et alj 120091 ). Whether CFs are present at 
very large radii (> 500 kpc) is currently unknown because 

of observational limitations. 

Cold fronts are not rare. IMarkevitch et al l ((2003) find 
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CFs in more than half of their cool core clusters, sug- 
gesting that most, if not al l such clusters harbor CFs 
(Markevitch fc Vikhlinii]|2007t ). Ichizzardi et al.l l|2010h sur- 
vey a sample of 42 clusters using XMM-Newton and find 
that 19 out of 32 nearby (2 < 0.075) clusters host at least 
one CF. Based on observational limitations, and orienta- 
tions of cold fronts, they conclude that this is a lower limit. 
While many CFs are directly attributed to mergers, out 
of 23 clusters in their sample that seem relaxed, 10 ex- 
hibit cold fronts. These 10 have systematically lower cen- 
tral entropy values. Another useful compilation of clusters 
for which co l d fron ts have been observed can be found in 
lOwers et al] l |2009h . Out of the 9 high-quality Chandra ob- 
servations of clusters with observed cold fronts, 3 of them, 
RXJ1720. 1+2638, MS1455. 0+2232 and Abel 2142 have non- 
disturbed morphologies. 

In what follows, we propose that some CFs are produced 
by merging of shocks. When two shocks propagating in the 
same direction merge, a CF is always expected to form. Its 
parameters can be calculated by solving the shock conditions 
and corresponding Riemann problem. Most of the scenarios 
in which shocks are produced (quasars, AGN jets, mergers) 
predict that shocks will be created at, or near, a cluster 
centre, and expand outwards (albeit not necessarily isotrop- 
ically). When a secondary shock trails a primary shock, it 
always propagates faster (supersonically with respect to the 
subsonic downstream flow of the primary shock), so colli- 
sions between outgoing shocks in a cluster are inevitable if 
they are generated within a sufficiently short time interval. 
We note that contact discontinuities can also form from head 
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on collisions between shocks, and various other geometrical 
configurations. We restrict ourselves in this paper to merg- 
ing of two shock propagating in roughly the same direction. 

In § [2] we solve for the parameters of CFs that are 
caused by merging two arbitrary planar shocks. We show 
results of shock-tube hydrodynamic simulation that agree 
with the calculated analytical prediction. In § U we de- 
scribe our spherical hydrodynamic simulations and realis- 
tic cluster profiles that are later used to investigate how 
these shocks may be produced, and discuss various ways to 
produce shocks in cluster settings. We identify two general 
mechanisms that create shocks in clusters: energetic explo- 
sions near the centre, or changes and recoils in the poten- 
tial well of the cluster that correspond to various merger 
events. These mechanisms are simulated in § 14.11 and § 14.21 
and also serve to demonstrate the robustness of the shock 
induced cold fronts (SICFs) that form. In §[5] we follow clus- 
ter growth and accretion over a Hubble time by using self- 
consistent ID hydrodynamic simulations of gas and dark- 
matter. The virial shock that naturally forms serves as a 
primary shock, and we demonstrate that SICFs form when 
secondary shocks propagating from the centre merge with 
it. We show that oscillations of this shock around its steady 
state position can produce a series of SICFs. In §[B]we dis- 
cuss aspects regarding the stability and sustainability of cold 
fronts in general, and derive some observational predictions 
that could serve to distinguish between CFs produced via 
this and other mechanisms. In § [7] we discuss our findings 
and conclude. Compact expressions for SICF parameters can 
be found in appendix § [X] 



2 MERGING OF TWO TRAILING SHOCKS 

In this section we examine two planar shocks propagat- 
ing in the same direction through a homogeneous medium 
and calculate the parameters of the contact discontinuity 
that forms when the second shock (secondary) overtakes the 
leading (primary) one. This problem is fully characterized 
(up to normalization) by the Mach numbers of the primary 
and secondary shocks, Mo and Mi. At the instant of col- 
lision, a discontinuity in velocity, density and pressure de- 
ve lops, corresponding to a Riemann problem (second case 
in lLandau fc Li fshitz 1959, §93). The discontinuity evolves 
into a (stronger) shock propagating in the initial direction 
and a reflected rarefaction wave, separated by a contact dis- 
continuity which we shall refer to as a shock induced CF 
(SICF). The density is higher (and the entropy lower) on 
the SICF side closer to the origin of the shocks. In spherical 
gravitational systems, this yields a Rayleigh- Taylor stable 
configuration if the shocks are expanding outwards. Now, 
we derive the discontinuity parameters. 



2.1 The Discontinuity Contrast 

We consider an ideal gas with an adiabatic index 7. Before 
the shocks collide, denote the unshocked region as zone 0, 
the region between the two shocks as zone 1, and the dou- 
bly shocked region as zone 2 (these definitions are demon- 
strated schematically in fig.[T]). After the merging, regions 
and 2 remain intact, but zone 1 vanishes and is replaced by 
two regions, zone 3 Q (the outer region for outgoing spherical 
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Figure 1. Schematics of shock merging. Pressure p is shown vs. 
an arbitrary spatial coordinate x, which is not fixed between the 
panels. Top: two shocks, propagating according to the arrows, in- 
duce transitions between zones -» 1 and between 1 — > 2. Middle 
panel, the instance of collisions between the shocks. The transition 
— > 2 does not correspond to consistent shock jump conditions. 
Bottom: at the Lagrangian location of the collision an isobaric 
discontinuity between 3 and 3; forms, separating between a for- 
ward moving shock (0 —¥ 3 ) and a reflected rarefaction (smooth 
transition between 2 and 3;). 



shocks; adjacent to zone 0) and zone 3i (inner; adjacent to 
2), separated by the SICF. We refer to the plasma density, 
pressure, velocity and speed of sound respectively as p, p, u 
and c. Velocities of shocks at the boundaries between zones 
are denoted by Vi, with i the upstream zone. We rescale all 
parameters by the unshocked parameters po and po- 
The velocity of the leading shock is 



= ito + Mo cq, 



where 



Pi 



1/2 
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for each zone i. Without loss of generality, we measure ve- 
locities with respect to zone 0, implying uo = 0. 

The state of zone 1 is related to zone by the Rankine- 
Hugoniot conditions, 
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Figure 2. Contrast q (labeled contours) of a CF generated by a 
collision between two trailing shocks with Mach numbers Mo and 
Mi, for 7 = 5/3. The color map shows the final Mach number 
Mf normalized by Mo Mi. 
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The equations are ordered such that each equations uses 
only known variables from previous equations. The state of 
the doubly shocked region (zone 2) is related to zone 1 by 
reapplying equations []J [5] with the subscripts 0, 1 replaced 
respectively by 1,2, and Mo replaced by Mi. 

Across the contact discontinuity that forms as the 
shocks collide (the SICF), the pressure and velocity are con- 
tinuous but the density, temperature and entropy are not; 
the CF contrast is defined as q = p3i/p3o- Regions and 3 
are related by the Rankine-Hugoniot jump conditions across 
the newly formed shock, 
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The adiabatic rarefaction from pressure P2 down to ps = 
P3i = P3o is determined by 
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The system can be solved by noting that the sum of eqs. (0) 
and ((8]) equals U2—U0, fixing p3 as all other parameters are 
known from eqs. |T]|5J). 

Finally, the Mach number of the new shock and the 
rarefacted density are given by 
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Closed-form expressions for q and Mf are presented in 
: [S] Figure [2] shows the discontinuity contrast q for vari- 
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Figure 3. The maximal SICF contrast (solid, left axes) as a 
function of 7. It occurs in a collision between a very strong shock 
Mo — > 00 and a trailing shock of finite Mach number Mi (dashed, 
right axes). For 7 = 5/3, q m ax = 2.65 with Mi = 6.65. 



ous values of Mo and Mi, for 7 = 5/3. Typically q ~ 2, 
ranging from 1.45 for Mo = Mi = 2, for example, to 
qmax = 2.653. The figure colorscale shows the dimension- 
less factor / = Mf /(M0M1), ranging between 0.75 and 1 
for 1 < {Mo, Mi} < 100. 

The maximal contrast q ma x depends only on the adia- 
batic index, as shown in fig. [3] It occurs for a very strong 
primary shock Mo —¥ 00 and a finite secondary Mi ; see § 
for details. For 7 = 5/3 we find g mal = 2.653, achieved for 
Mo > 1 and Mi ~ 6.65. 



2.2 Planar Hydrodynamic Example of SICF 
Formation 

To demonstrate the formation of an SICF we present in fig.|4] 
results from a ID planar hydrodynamic code. The numeri- 
cal scheme of the hydrodynamic calculation is Lagrangian, 
with the locations and velocities defined at edges of cells 
and thermodynamic variables defined at the centres of cells. 
The density p and internal energy e, are the free thermody- 
namic variables, setting the pressure p, temperature T and 
the entropy S according to an ideal gas equation of state. 
Throughout the paper we consider monoatomic gas, with 
an adiabatic constant of 7 = 5/3 unless stated otherwise. 
The time propagation is calculated by a leap-frog explicit 
scheme. The timesteps here are adaptive, and set by a stan- 
dard Courant Friedrichs Levy condition. Discontinuities are 
integrated over by employing first and second order Von 
Neumann artificial viscosity. 

The simulation presented in fig.|4]was run with 500 cells. 
Without loss of generality, we can assume po = eo = To = 1. 
The definition of temperature here implies an arbitrary heat 
capacity that does not enter into the equations of motions. 
We define the entropy as S = T/p 7-1 (the exponent of the 
usual formal definition) so So = 1 as well. We use external 
boundary conditions of a piston propagating into the mate- 
rial from X — in the positive direction at Mach number 
Mi = 2.35. At time t = 0.04 the piston's velocity is increased 
to create a second shock propagating in the same direction, 
with M2 = 1.76 with respect to the previously shocked gas. 
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Figure 4. Shock tube simulation of a contact discontinuity 
(SICF) formed by merging of two shocks. The piston propagates 
from x = in the positive direction at a supersonic speed cre- 
ating the shock between — > 1. At t=.04 the piston's velocity 
abruptly increases causing a second shock to propagate between 
1 — ¥ 2. At the position and time where the shocks merge a con- 
tact discontinuity forms between regions 3 and 3;. Top: density 
(normalized so that po = !)■ Red lines mark the trajectories of 
some Lagrangian cell boundaries. The zones here are analougous 
to those in § [5] and fig. [TJ Bottom: entropy (normalized so that 
So = To/pg = 1). Note that there is no entropy change between 
zone 2 and 3i because the gas rarefies adiabatically between them. 
The contrast of the SICF, q, is the density jump between zones 
3 and 3i. 

The contact discontinuity that forms has a density ratio3 
q = p3i/p3o = 1.466 while the theoretical prediction accord- 
ing to equationsfTlllOlis q = 1.459 - in satisfactory agreement 
for the resolution used, and for the accuracy needed in the 
discussion presented here. The same technique has also been 
performed for larger Mach numbers and yielded results with 
comparable accuracy. 



The SICF strength is calculated in the simulation according to 
q = (Szi/ Sso) -1 / 1 because in the rarefied zone, 3i the density is 
not constant, and so harder to measure directly, while the entropy 
is constant throughout the zone. 



3 NUMERICAL INVESTIGATION OF 
SPHERICAL COLD FRONTS 



3.1 Driving Shocks 

We explore different ways to drive central shocks through 
clusters of galaxies. The first, most intuitive method is ex- 
ternal energy injection at the centre (i.e. explosions), which 
would correspond, for example, to mechanical energy re- 
leased in AGN bursts. Observations in X-ray a nd radio indi- 
cate a wide variety of sho cks and sound waves (|Fabian et al.l 
120031 ; iForman et al.1 2007) emitted repeatedly, perhaps even 
periodically from a central AGN. Many clusters also include 
hot and diffuse bubbles, that are most likely inflated by 
AGN jets that a re decollimated by the amb ient cluster gas 
(see review by iMcNamara fc Nu lscn 2007, and reference 
therein). If these bubbles are a plausible mechanism to coun- 
teract the overcooling problem, the energies associated with 
them must be compar able to the cooling rate of the clus- 
ters (> 10 44 erg s ec" 1 ; lEdge et al.|[l990l ; TDavid et al J 1 19931 ; 
lMarkevitchlll998T ) over some fraction of the Hubble time, 
indicating that the total energy injection is > 10 61 erg. The 
process of inflation of these bubbles sends shocks through 
the ICM, but at the same time the extremely hot, under- 
dense bubbles alter the radial structure of the ICM in a non 
trivial manner. In § 14.11 we demonstrate (noting the limi- 
tations of the spherically symmetric analysis in that case) 
that mergers between shocks driven by such central explo- 
sions can produce SICFs. 

Shocks can also be driven gravitationally. Mergers of 
galaxies and subhaloes change the structure and depth of the 
potential well causing sound waves and shocks to propagate 
through the haloes. These shocks typically do not consid- 
erably disturb the ICM structure, leaving a static gaseous 
halo after they propagate. We show in § 14.21 that abrupt 
changes to the gravitational potential of the halo directly 
produce shocks. In particular, in the first passage stage of a 
merger, two outward propagating shocks are produced: one 
when the halo contracts due to the increase of the gravi- 
tational force, and another when the halo abruptly rarefies 
once the merging subhalo flies out. The numerical investiga- 
tion of explosion and merger driven SICF is performed us- 
ing a spherical one dimensional hydrodynamic code within 
a static pote ntial well (a version of t he code "Hydra" doc- 
umented in iBirnb onn fc Dekellliooi that is stripped from 
the dynamic dark matter evolution and self-gravitating gas) 
that, along with the initial conditions, is described in § 13.21 

In §[5]we investigate the interaction of a secondary shock 
with the virial shock that is always expected at edges of 
clusters. To produce this primary shock, we simulate the 
cosmological evolution of a cluster by using a hydrodynamic 
scheme that includes a self-gravitating gas that flows within 
self-gravitating, dynamic dark matter shells, and starts from 
spherical cosmological overdensities at a redshift z = 100. 
We show that perturbations can either be invoked manually, 
or occur naturally when infalling dark matter has sufficiently 
radial orbits so it interacts with the central core directly. In 
this context, we also show that low amplitude perturbations 
in the core send sound waves that can steepen into a shock 
that interacts with the virial shock. 
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3.2 The "Hydra" Hydrodynamical Code 

The cosmological implications of shock interactions in clus- 
ters are tested using ID, spheri cal hydrodynamic simula - 
tions performed with "Hydra" |Birnboim fe Dekell 1200 j) . 
The code is run here in two modes: cosmological and static 
DM (dark matter) modes. In §[4]we use the code in the static 
DM mode, a simplified configuration in which the baryons do 
not self-gravitate, and there is no evolution of the dark mat- 
ter. Rather, baryons start in hydrostatic equilibrium within 
a fixed cluster-like potential well. In §[5]we use the code in its 
cosmological mode, utilizing the full capabilities of the code 
to evolve baryons and dark matter self-consistently from a 
cosmological initial perturbation at high redshift to z — 0. 
Next, we summarize the physical processes implemented in 
both modes of "Hydra" and then highlight the difference in 
boundary conditions and initial conditions corresponding to 
each modus operandi. 

Similar to the planar code described in § 12.21 the code 
uses a Lagrangian scheme. The positions of baryonic spher- 
ical shells are defined by their boundaries (with the inner- 
most boundary at r = 0), and the thermodynamic properties 
are defined at centres of shells. A summary of the fluid equa- 
tions that are solved, and the numerica l difference scheme 
is available in iBirnboim fc Dekell (|2003h . The code (in its 
cosmological mode) also evolves thin discrete dark matter 
shells, that propagate through the baryons, and interact 
with them gravitationally. The full discrete force equation 
for the "i"th baryonic shell is thus: 

2 AP t GMj 1 
U - ^ AMi {n + eY + ( ' 

with ri, Mi,ji the radius of boundary i, the mass (baryonic 
+ DM) enclosed within that radius, and the specific an- 
gular momentum prescribed to this shell respectively, and 
APi, AMi the pressure difference and baryonic mass differ- 
ence between the centres of shells i + 1 and i. e is the grav- 
itational smoothing length (50 pc and 500 pc in the static 
DM and cosmological modes respectively). 

The dark matter and baryonic shells are assigned an- 
gular momentum that acts as an outward centrifugal force 
designed to repel shells that are close to the singularity at 
the centre, in analogy to the angular momentum of a single 
star, or of the gaseous or stellar disk in realistic systems, 
keeping them from falling into the centre of the potential 
well in the absence of pressure support. In the cosmologi- 
cal mode the shells are initially expanding with a modified 
Hubble expansion, and the angular momentum is added as 
each shell turns around, corresponding to the physical gen- 
eration of angular momentum via tidal torquing around the 
maximal expansion radius. In the static DM mode, the an- 
gular momentum is present from the beginning of the sim- 
ulation, and is taken into account when the hydrostatic ini- 
tial profile is calculated. The angular momentum prescrip- 
tion of the baryons is determined by the requirement that 
as gas passes through the virial radi us, its angular mom en- 
tum will be some fraction (A = 0.05, iBullock et afl|200ll ) of 
the value needed to support it on circular motions. As ex- 
pected, the gas only becomes angular momentum supported 
near the centre, and when gas cooling is turned on (in the 
simulations presented below, the angular momentum of the 
baryons is never important, and is presented here only for 



completeness). The dark matter is prescribed angular mo- 
mentum according to the eccentricity of a shell's trajectory 
as it falls in through the virial radius: j = v v j T r vlr \/2(l — /3), 
(corr esponding to: Vg — «? = (1 — /3)t^; [Bmnev fc Tremaing 
1 19871 . § 4.2). w v i r , r v i r are the virial velocity and radius, and 
v r ,vg,v,p are the radial and two tangential velocity compo- 
nents of an infalling trajectory. Low values of /3 indicate 
large angular momentum, and /3 = 1 corresponds to purely 
radial orbits. j3 is constant for all the dark matter shells in 
the simulation. Gas cooling ca n be turned on, by interpolat - 
ing from a tabulated version of lSutherland fc Dopital {l993), 
with a prescribed metallicity. 

Rather than use a leap-frog integration scheme, which 
is sufficient for most explicit hydrodynamic schemes, we use 
here 4th order Runge-Kutta integration over the coupled 
dark matter and baryonic timestep. This allows us to derive 
a timestep criterion for the dark matter shells, by compar- 
ing changes in the velocities and radii of the dark matter 
and baryons between the 1st and 4th order propagation 
and limit this difference to some fraction of the velocity 
and radius. The timestep is determined by the combined 
Courant Friedrichs Levy condition and this requirement. If 
the criterion is not met, the code reverts to the beginning of 
the timestep and recalculates the next step with a reduced 
timestep. 

The numerical integr ation scheme, angular mom entum 
and cooling is described in IBirnboim fc Dekell l|2003h . In ad- 
dition, we report here on two new ingredients to the cal- 
culation: adaptive changes to the grid, and ID convective 
model. Baryonic shells are now split and merged if they be- 
come too thick or thin respectively. A shell is split if its 
width is larger than some fraction (typically 0.2) of its ra- 
dius, or if it is larger than some fraction (typically 4) of 
the shell directly below it and above it, provided that it is 
larger than a minimal value (typically 2 kpc), and that the 
shell is not in the angular momentum supported region near 
the centre of the halo. Shells are merged when the cumula- 
tive width of two adjacent shells is smaller than some value 
(typically 0.1 kpc). Shell merging is helpful in cases where 
two shells are pressed together, decreasing the timestep to 
unreasonable values (reasonable timesteps for cosmological 
simulations are above 10 -7 Gyr). An additional component 
included in "Hydra" for the first time is a ID mixing length 
convection model, that acts to transfer energy out wards in 
regio ns where the entropy profile is non- monotonic jSpiegel 
1963). This model is described in detail in lBirnboim fc Dekel 
l|2010r) . and is used here only in the final case described in 
§ [5] in its maximal form: bubbles can rise until they reach 
the local speed of sound. In this work, the inclusion of the 
maximally efficient mixing length theory convection ensures 
that the entropy is always monotonically increasing, in an 
energetically self-consistent way. 

The initial conditions for a cosmological run are derived 
by requiring an aver age mass evolution h istory for a halo 
with some final mass (Neisteinetal ] |2006T ). The overdensity 
of each shell at the initial time is calculated so that the shell 
would pass through its virial (riso) radius at the required 
time (the procedure is described in detail in lBirnboim et ail 
l2007f ). This procedure has an additional benefit of setting 
the initial grid in a way that would ensure that shells will 
pass through the virial radius at constant time intervals. In 
this procedure, and in the cosmological simulation, a force 
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Figure 5. Radial profiles of cluster initial conditions. Top: den- 
sity (red, left axis) and total enclosed mass (green, right axis) 
of the cluster. Bottom: Temperature (red, left axis) and entropy 
(green, right axis), as a function of radius (upper x axis) and virial 
fraction (lower x axis). The total mass profile, and the tempera- 
ture profile are prescribed, as well as a pivot point in the entropy 
(value of 10 keV cm 2 at r = 0.01i? v ir)- The density and entropy 
profiles are calculated so that the gas is in hydrostatic equilib- 
rium. The temperature profile is set to exhibit a decrease by a 
factor ss 1.5 between 0.1 to 0.01-f? v i r . 



corresponding to the cosmological constant has been added 
to eq. 111 ! making the simulat ion completely consistent with 
ACDM ijBirnboim et al.ll2007l ) (the cosmological parameters 
used throughout the paper are: fi m — 0.3, £Ia = 0.7, h = 
0.7). This has little overall effect, and essentially no impact 
once the shells are collapsed, and the overdensities become 
non-linear. 

In the stripped down, static DM mode, there are no 
dark matter shells, and the gas does not self-gravitate. In- 
stead, the value M; of eq. [TT] is interpolated from a prede- 
fined lookup table. In this mode, a rigid external boundary 
condition is applied, in contrast to a vacuum boundary con- 
dition in the cosmological mode. The inner boundary con- 
dition in both cases is trivially satisfied because all terms in 
eq. [11] are individually zero there. 

For the static DM simulations we wish to define initial 
conditions which are typical of those in clus ters. The total 
mass profile is defined by an NFW profile (jNavarro et al.l 
119971 ). with M v i r = 3 x 10 14 M Q , and conce ntration set 
to c = 5.4 according to iBullock et all (|200lT ). The virial 
radius is 1730 kpc. A central ga l axy (BCG) is modeled 
by a Hernquist profile (|Hernquistl Il990l ) superimposed on 
the NFW profile, with a sharp cutoff at 10 kpc and 
concentration a = 6 kpc normalized to Mbog — 3 x 
1Q 11 Mq. We define an inward decreasing temperature pro- 
file i|Leccardi fc Molendill200Sl ) that peaks at tt = 0.1 R vil 
at T(rr) = Tvir = 2.7 x 10 7 K and drops inwards according 
to T(r) = T v i r (r /rr) ' 2 ■ Outwards, it is linearly decreas- 
ing from T(r T ) = T vit to T(R vil ) = 0.7 T vir . Once the en- 
tropy of a single point is defined, the density profile can be 



constructed by the hydrostatic requirement without further 
degen eracy. We set S(0.01R v i T ) = 10 keV cm 2 as our pivot 
point (|Cavagnolo et al1l2009l ). The final baryonic mass of the 
resulting profile in the ICM is 2.5 x 1O 13 M , which is w 8% 
of the total mass. The density, total mass, temperature and 
entropy profiles are plotted in fig. [S] We note that while the 
constructed initial profile reasonably resembles that of a re- 
alistic cluster, the mechanism that is studied here, of SICFs, 
is not sensitive to these details; any merging between shocks, 
in any profile, should produce them. 



4 STATIC DM SIMULATIONS 
4.1 Cluster Explosions 

The first mechanism we investigate to produce shocks is that 
of cosmic explosions. Many processes can cause abrupt en- 
ergy injection at centres of clusters, among which are AGNs 
or starbursts at the BCG and galaxy-galaxy mergers. The 
discussion in this section is confined to injection to the bary- 
onic component itself (changes to the gravitational potential 
well will be discussed next). As shocks sweep through mat- 
ter, they constantly loose energy to heating the gas, and 
when the shocks expand spherically, they also weaken due 
to the geometrical increase in the surface of the shock 0. We 
seek here to create shocks that will survive to some consid- 
erable fraction of the halo radius, requiring that the amount 
of energy that is injected will be comparable to the total 
thermal energy of the gas in the cluster. This energy scale 
is similar to that required to solve the overcooling problem, 
indicating that any solution to the cluster overcooling prob- 
lem via strong shocks will produce shocks that will expand 
to a significant fraction of the halo radius. 

Fig. [6] describes the time dependent evolution of an ini- 
tially hydrostatic cluster atmosphere, as two explosions are 
driven through it. The simulation has 1, 000 logarithmically 
spaced shells, between 0.0005 R v i v and R V1I . Here, and in the 
following examples, the cooling is turned off, which is justi- 
fied in the absence of any feedback mechanism in the simula- 
tion that would counteract the overcooling. The explosions 
are driven by setting a homologic (v — v axp r /r ex p) veloc- 
ity profile with highly a supersonic velocity (v cxp ), sharply 
cutting off beyond some radius r exp . For the first explo- 
sion, these parameters are v exp = 10,000 km sec" 1 and 
r exp = 20 kpc, injecting 1.7 x 10 61 erg into the core of the 
cluster. The second explosion, 3 x 10 r yr later, has a veloc- 
ity of u CX p = 20,000 km sec -1 at r exp = 40 kpc, and in- 
jects 4.7 x 10 61 erg. These large values are not a-typical for 
the energy estimated in ra dio bubbles l|Brrzan et alj [20041 ; 
iMcNamara fc Nulsenl [20071 ). While no buoyancy is possible 
within a spherical symmetric framework, the explosions do 
create extremely hot and dilute regions, and drive strong 
shocks. The two shocks can be traced in fig.|S]by the breaks 
in the plotted Lagrangian lines (velocity jump) , by the tem- 
perature jump, and by the regions in the simulation with 
large artificial viscosity (black dots). At the location of 

2 IWaxman fc Shvartj dl993r i; iKushnir etall j2005h showed that 
for power law density profiles and strong shocks, as long as the 
mass diverges with radius, or, p gas oc r~ m , with w ^ 3, the shock 
will weaken 
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Figure 6. Time evolution of an initially hydrostatic halo show- 
ing explosion driven shocks from the centre propagating outwards. 
The energy injection parameters are described in the text. Top: 
entropy color map on the radius-time plane. The red thin lines 
mark the position of every 50th Lagrangian shell boundary. The 
black dots mark shocks, defined as artificial viscosity pressure 
exceeding 10~ 11 erg cm~ i . These shocks are driven by two explo- 
sions, at t = 0.1 Gyr and 0.f3 Gyr. The entropy of the shocked 
gas is saturated in this colormap, but the entropy jump of the 
SICF is visible at the Lagrangian location of collision between 
the shocks, middle: temperature colormap of the same simula- 
tion. The temperature increase induced by the shocks, and the 
temperature jump at the SICF (by a factor of f .3) is visible. Bot- 
tom: Pressure colormap. The pressure across the SICF is smooth, 
as expected for cold fronts. 



merging between the shocks (at t = 0.15 Gyr and radius 
~ 150 kpc), an SICF (contact discontinuity) forms (visible 
in the entropy colormap) that propagates with the matter. 
We have followed the simulation for 1 Gyr and, as expected 
in the absence of any diffusive mechanisms, the SICF retains 
its strength. 

The energies and time interval between the shocks 
have been chosen arbitrarily to create a shock collision at 
~ 150 kpc. A series of sharp, frequent and strong cen- 
tral shocks is expected in some self- consistent radiative and 
mechanical AGN feedback models (|Ciotti fc Ostrikerl f2007; 
ICiotti et al.l 120091 ). and is indicated in some observations 



(for example, iFabian et all I200&1 ; iForman et all 120051 . for 
Perseus and M87 respectively). Non-spherical, frequent col- 
lisions may occur between the shocks produced by the two 
bubbles of a radio active AGN, at a position roughly above 
one of the bubbles (the generation of two shocks correspond- 
ing to observed bubbles, as well a s radial soundwave s that 
steepen into a shock is discussed in IFabian et aill2003l ). The 
shock from the near-by bubble (primary) arrives to that po- 
sition first, and the shock from the opposite bubble (sec- 
ondary) arrives later. Given our ID framework, we do not 
attempt to map the parameter space of possible explosions. 
Since shocks weaken as they propagate outwards, weaker 
shocks which occur within a shorter time interval will pro- 
duce smaller radii SICF with the same strength as stronger 
shocks driven within a larger time interval. 

4.2 Gravitational Perturbations: Merger Induced 
Shocks 

If the core of an otherwise hydrostatic cluster is gravitation- 
ally perturbed, the reaction of the gas to the perturbation 
often results in shocks. While explosions, discussed in § 14.11 
typically considerably heat the shocked gas near the cen- 
tre, gravitational perturbations make a more subtle change 
to the gas, while forming strong outward expanding shocks. 
The core is perturbed when gas or dark matter is accreted 
to the centre of the cluster, either smoothly (creating weak 
perturbations) or by mergers, (causing a more violent and 
abrupt shock). In this section we describe the shocks that 
occur as a result of a large change to the core potential depth 
of a cluster due to mergers. In §[5] we show that weak oscil- 
lations in the core send sound waves through the gas, that 
steepen into shocks. 

Fig. [7] shows a merger event of a cluster with the same 
initial conditions described in § 13.21 and fig. [5] with a dense 
substructure with total mass of 3 x 1O 13 M (a 10:1 merger 
ratio). The simulation describes, in a simplified form, how 
the main halo will react to the first passage of the substruc- 
ture near its core. During its first pericenter passage, the 
subcluster will first cause the core to contract, and then, as 
it flies out, allow the core to relax again. Here, we model 
this in a simplistic way by adding, and then removing (after 
2 x 10 7 yr), a massive object at the centre of the simulated 
cluster. The addition of the mass causes all the cluster's at- 
mosphere to contract, affecting the inner region more than 
the outer. A coherent flow inwards begins, that stops as an 
outward expanding shock passes through the cluster, read- 
justing it to a new hydrostatic equilibrium (note in fig.[7]that 
after the passage of the primary shock the post-shock gas is 
almost at rest). Once the central object is removed, the halo 
jolts out from the inside out (again, because the addition 
of the central force affects the inner core the most), caus- 
ing another shock to occur. The energy for these two shocks 
is taken from the orbital energy of the merging subcluster, 
and would cause it to spiral in. In the figure both shocks 
are clearly visible in the temperature maps, and through the 
black dots corresponding to strong artificial viscosity. At the 
instant of merging between the shocks (at t = 0.044 Gyr and 
radius ~ 25 kpc), an SICF forms which is seen in the en- 
tropy and temperature colormaps. The time interval for the 
merger interaction is chosen manually, but does not seem 
improbable noting that an object approaching the centre 
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Figure 7. Same as fig. [6] but now the potential well is modified 
by artificially adding a central object of mass 3 X 10 13 Mq, which 
is then removed 2 X 10 7 yr later (corresponding to a first passage 
of a 1:10 ratio subhalo). Note that the scales are different from 
fig- 111 

at 1000 km sec -1 , roughly the virial velocity, passes 20 kpc 
during that period. Shocks associated with sequential core- 
passages of an inspiraling subhalo could also induce large 
scale SICFs. While shoc ks produced by merger events have 
been seen in simulations (|McCarthv et al]|2007l ). the forma- 
tion of SICF in these simulations hasn't been tested yet, and 
is left for future work. 



5 VIRIAL SHOCKS AND SICFS 

At the edges of galaxy clusters there are virial shocks with 
typical Mach numbers ~ 30 — 100, heating the gas from 
~ 10 4 K to ~ 10 7 K by converting kinetic energy into ther- 
mal energy. Virial shocks are probably ob served by Suzaku 
|George et all 120091 : Irloshino et"aHl2010l ). The rate of ex- 
pansion of a virial shock is set on average by the mass 
flux and velocity o f the infalling material (see for example 
iBertschingerl 1 1985(1 Secondary shocks that originate from 
the centre of clusters, due, for example, to the mechanisms 
discussed in §[3] will collide and merge with the virial shock, 
creating SICFs at the locations of the virial shock at the 
corresponding times. Since shock mergers with virial shocks 
require only one additional shock to form (rather than two 
in the general cases described in § H), and no particular 



timing is required, such SICFs are perhaps even more prob- 
able. When the Mach number of the primary (virial) shock 
is Mo > 3, which is always the case, the strength of the 
SICF becomes almost independent of the secondary shock, 
allowing for a rather narrow range of strengths around q ~ 2 
as long as the secondary shock is sufficiently strong (fig. 0. 

We simulate the merging of a virial shock with a sec- 
ondary by evolving self consistently the dark matter and 
baryons from an initial perturbation in the cosmological 
mode described in § 13.21 These initial conditions provide 
a reasonable laboratory for gas dynamics in the ICM, pro- 
ducing correct temperatures and densities, and tracking the 
virialization process of the gas in a self consistent way. In 
order to capture the production of a secondary shock and 
its interaction with the virial shock in ID simulation, we 
examine various artificial schemes which mimic the full 3D 
dynamics, as discussed below. A more detailed simulation 
of cluster evolution would require 3D simulations that in- 
clude AGN feedback and merger and smooth accretion for 
baryons and dark matter. The cosmological simulations here 
have been performed with 2, 000 baryonic shells, and 10, 000 
dark matter shells, sufficient according to our convergence 
tests. 

For the average mass accretion histories used here 
(§ 13. 2p . the accretion rate grows with halo mass, even af- 
ter ta king into account the late time decline in accretion 
rates (|Dekel et al.l[20 09), so shells that fall later are more 
massive. When late massive dark matter shells fall into the 
inner core, their mass becomes comparable to that of the 
gas enclosed within that shell, causing it to vibrate stochas- 
tically. These core vibrations drive sound waves that are 
emitted from the centre and steepen to create shocks (the 
energy for these shock is taken from the orbits of the dark 
matter shells) . We can control the level of stochastic noise by 
increasing or decreasing the angular momentum prescribed 
to the dark matter. In the first two "synthetic" examples 
shown here, we use /3 — 0.7 (see § I3.2|l to reduce noise in 
the evolution, and then manually add perturbations at time 
— 7 Gyr (z = 0.84). In the last example, we decrease the 
dark matter angular momentum (/? = 0.991), and use the 
numerical stochastic noise as a proxy for the perturbations 
that halo cores undergo during their formation and accre- 
tion. We demonstrate that this noise naturally give rise to 
shocks that interact with the virial shock. This interaction 
seems to be quasi periodic, and we suggest a mechanism 
for driving this behavior. In this final example, the entropy 
profile is sometimes non-monotonic, implying that convec- 
tion might be important. We test the effect of convection 
using our ID mixing length theory and show that it does 
not significantly reduce the strength of the SICFs that are 
produced, and does not change the overall evolution signifi- 
cantly. 

5.1 Manual Initiation of Secondary Shocks 

Fig. [8] describes a typical evolution of a cluster that ends up 
with mass 3 x 1O 14 M0 at z = 0, and subsequent interaction 
of the virial shock with a secondary merger induced shock. 
The merger here is 10:1, with a central mass of 3 x 1O 13 M0 
instantaneously added at —7 Gyr (at that time the cluster's 
virial mass is ~ 10 14 Mq). Unlike the procedure described in 
§ 14,21 we do not remove the merged mass, creating only one 
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Figure 8. Cosmological mode simulation of the evolution of a M(z = 0) = 3 X 10 14 Mq cluster. The gray lines mark the evolution of 
every 50th Lagrangian shell, and the colormaps are entropy (top) and temperature (bottom). The initial Hubble expansion, turnaround 
and the virial shock are visible. Initially, at z = 100, shells expand due to the near-Hubble flow. At time —7 Gyr a central object with 
mass 3 X 1O 13 M0 is permanently added, causing a shock to propagate outwards. The collision between this shock and the virial shock 
creates an SICF visible in the entropy colormap which persists to z = 0. 



shock. The Hubble expansion, turnaround and virial shocks 
are seen, as well as the SICF formation at time ~ —5.5 Gyr 
and at radius « 700 kpc. After the deepening of the grav- 
itational potential well, the secondary shock mediates the 
re-establishment of hydrostatic equilibrium. Thus, the post- 
shock gas, and the SICF, are roughly at rest in the Eulerian 
laboratory frame of reference. The post-shock temperatures, 
particularly below 100 kpc are unrealistically high, corre- 
sponding, perhaps, to disturbed morp hologies seen in clus - 
ter merger events (the bullet cluster; iTucker et al.|[l998l ). 
although no entropy inversion occurs there. This is remedied 
somewhat when the shock is initiated by a series of weaker 
shocks and sound waves that steepen into a strong shock 
further from the cluster's centre (fig. , and is completely 
alleviated when the shocks are due to persistent stochastic 
noise near the core (§E2J, generating profiles which are also 
consistent with quiescent cool core clusters. 

We produce a weaker, more local perturbation by en- 
gineering a perturber that sends weak shocks and sound 
waves that steepen into a shock. Fig. [9] shows an evo- 
lution of the same cluster as in fig. [H] but with a cen- 
tral perturber who's mass is sinusoidally changing between 
-7 < t < -6.75 Gyr with an amplitude of 3 x 10 13 M Q and 
a period of of 5 x 10 7 Gyr. After that, a constant mass of 
5 x 10 12 Mq is left. The steepening of the shocks is visible, 
and the temperature increase of the gas by the secondary 



is much smaller, because, contrary to the permanent mass 
addition in fig. [HI no net contraction is caused. This example 
serves to show how shocks can form from "noise" generated 
near the core, but still heat the centre from the inside out, in 
a way that would correspond to a merger event but not to a 
quiescent cool core cluster. The entropy profile here becomes 
non-monotonic near the centre. 



5.2 Self Consistent Initiation of Secondary Shocks 

We allow numerical noise to perturb our hydrodynamic sim- 
ulation by reducing the angular momentum of dark matter 
so that late infalling, massive dark matter shells penetrate 
the cluster core. The resulting time evolution and final pro- 
file are shown in figures \W\ - 1121 As soon as the cluster 
progenitor forms, at high redshift, sound waves and weak 
shocks begin propagating through the cluster, as the gas 
constantly adjusts to the varying potential well near the 
centre. The weak shocks steepen and accumulate to create 
strong shocks that propagate outwards and merge with the 
virial shock, creating an SICF on each merger instance. From 
fig. [10] it is evident that the shocks are weak, and that they 
do not alter the inner profile considerably. 

The shocks that are invoked by the stochastic noise 
seem to be periodic, with a period corresponding to about 
twice the sound (or shock) crossing time through the cluster. 
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Figure 9. A cosmological mode simulation with the same initial conditions as in fig. [8] but with the core perturbed briefly between 
—7 and —6.75 Gyr (see § 15.11 for details). The perturbation is a sinusoid in time, with an amplitude of 3 X 1O 13 M0, and a period of 
5 X 10 7 Gyr, following by a constant left-over mass of 5 X 1O 12 M0. Weak shocks and sound waves are sent outwards, steepening into a 
single shock that collides with the virial shock, producing an SICF. 



Lacking a formal analysis for this oscillatory mode, we sug- 
gest, based on these simulations, the following explanation. 
The location of the virial shock is regulated such that the 
pressure below (towards the cluster's centre) wi ll support the 
shock and the ram pressure of the infalling gas (|Bertschingeil 
1 19851 ). While this process would imply a smooth expansion 
of the shock, the virial shock can oscillate around its steady 
state trajectory in response to perturbations in the halo. 
When the post-shock gas is over-pressurized, the virial shock 
adjusts itself by expanding faster, relaxing this overpressure. 
If the relaxation overshoots, the pressure would drop below 
its steady state value, causing the shock to decelerate. The 
halo thus oscillates between these two phases depending on 
the virial shock's velocity: in one it is overly slow, produc- 
ing pressurization of the post-shock gas, and in the other 
it is overly fast, causing rarefaction. These modes are inter- 
laced. A rarefaction cycle begins with the secondary hitting 
the virial shock sending rarefaction waves inwards that are 
reflected at the centre and propagate outwards, ultimately 
reaching the virial shock making it decelerate. A compres- 
sion, or secondary shock cycle collects compression waves 
that are transmitted inwards when the virial shock is too 
slow, and are reflected through the centre, ultimately steep- 
ening into what would become the secondary shock. Circum- 
stantial evidence for this mode can be seen in fig. 1111 where 
the shocks are automatically traced, and the trajectory of 
one rarefaction wave reflected from a secondary- virial shock 



merger is highlighted by hand (dashed curve, based on a 
high resolution time sequence). Note that this periodicity 
is also seen (albeit much more weakly) in the two manual 
excitations of a secondary shock described in figures [8] and 
[9] The periodic SICFs that are produced can be seen in the 
entropy profiles in the upper panel of fig. 1101 in the lines 
marking the gradient of the entropy in fig. 111! and in the fi- 
nal profile of this simulation shown in fig. [12] The shocks and 
rarefactions propagating though the halo, and the alternat- 
ing speed of the viri al shocljfl are observed in th e 3D galactic 
halo simulations of [Keres fc Hernquist! (|2009l rl. Note that 
this qualitative argument is similar to the one described in 
§ [2] The overpressure there is caused by a secondary shock 
that perturbs the equilibrium solution of a single shock: the 
primary accelerates, and a rarefaction is sent back, relaxing 
the gas. 

The calculations presented in Figures 1 101 121 are adia- 
batic. When cooling is turned on, and in the absence of 



3 A movie of the evolution of radial profiles in time is available 
at http : //www . cf a. harvard. edu/~ybirnboi/SICF/sicf .html 

4 Dusan Keres, (private communication). Propagation of 
merger induced shocks and the subsequent bounce of the 

merge 



virial shock as these secondaries 



with it are seen 



http : //www. cf a. harvard. edu/~dkeres/movies/Blhr_10nl28_large_ 

There, left panel is gas density, right panel temperature, and a 
particularly notable shock forms at z 1.7. 
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Figure 10. A cosmological mode adiabatic simulation of cluster evolution with initial conditions as described in the § 15,21 The angular 
momentum of the dark matter is sufficiently small to allow massive, late accretion of dark matter shells to propagate and reach the 
centre, increasing the amount of stochastic "noise" with respect to the levels seen in figures \E\ and [9] A series of secondary shocks merge 
with the virial shock, creating a series of SICFs. 



feedback, this simulation suffers from overcooling, creating 
an overmassive BCG of 2 x 1O 12 M0 and BCG accretion rates 
(corresponding to star formation rates) of > 1 00 Mpyr -1 , 
and the luminosity exce e ds the L x —T rela tion (|Edge et al.l 
Il990l ; iDavid et al][l993l ; lMarkevitcblll99ct ). The excitation 
of the basic mode, and periodic CFs, are observed in all 
these simulations. Since the virial shock's strength is non- 
monotonic, it is not surprising that between SICF formation, 
the post-shock entropy is non-monotonic. This introduces a 
problem in ID simulations where convection due to entropy 
inversion cannot naturally occur. We test this by rerunning 
the adiabatic simulat ion with ID convection according to 
mixing length theory (jSpiegelll 19631 ). with the mixing length 
coefficient set such that the convection is limited only by 
the requirement that bubbles never exceed the sonic speed. 
This acts to redistribute the energy and entropy between 
SICFs, but does not change the overall dynamics. The final 
profile in this case is plotted in the dashed lines of fig. 1121 
The strength of the SICF in this case is reduced, but it is 
still clearly visible, with typical values of q « 1.5. It is highly 
improbable that convection operates in its maximal physi- 
cal eff iciency, particularly i n the presence of ICM magnetic 
fields. IParrish et all (|2008l . 120091 ) show that in the presence 
of weak magnetic fields the effective convection is smaller by 
many orders of magnitude than the heat fluxes carried by 
conduction, implying that it is highly suppressed from its 



mixing length theoretical value (Ian Parrish, private com- 
munication) . 

We note a qualitative similarity between reverberations 
of the virial accretion shocks discussed here, and another 
prominent case of accretion shocks in astrophysic s - of type 
II core collapse supernovae. iBurrows et all (2007) find that 
stalled accretion shocks around type II cores are unstable to 
2D (g-mode) reverberation that, after enough cycles, accu- 
mulate sufficient energy and amplitude to cause an outbreak 
of the stalled shock. While the standing accretion shock in- 
stability (SASI) is predominantly 2D, and acts on differ- 
ent scales than discussed here, it does draw its energy from 
the accreted gas, and is perturbed by sound waves that are 
emitted from the vibrating core, that accumulate to a single 
unstable mode. 



6 STABILITY AND OBSERVABILITY OF 
SHOCK INDUCED COLD FRONTS 

6.1 Stability 

The stability of CFs limits the time duration over which 
they are detectable, and so is important when comparing 
CF formation models with observations. Various processes 
can cause a gradual breakup or smearing of the CF. They 
act in other CF formation models as well. 

Thermal conduction and diffusion of particles across 
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Figure 11. Evolution of a galaxy cluster from cosmological ini- 
tial perturbations. The final mass of the cluster at z = is 
3 X 1O 14 M0. The radius and time of Lagrangian shells (every 
25 th shell) are plotted in red thin lines. Shocks are traced by 
large Lagrangian derivatives of the entropy (dlnS/dt > .1 Gyr — 1 , 
blue dots), and CFs are traced by their large entropy gradients 
(dlnS/dlnr > 0.5, green dots). Rarefaction waves can be seen as 
small motions of the Lagrangian shells (illustrated by the dashed 
magenta curve, manually added based on a time series analysis). 
Their trajectory is approximately the reflection of the preceding 
outgoing compression, time inverted about the last secondary- 
virial shock collision. 




Radius [kpc] 

Figure 12. Thermodynamic profile of the simulated cluster 
shown in figures 1 1 01 1 1 1 at z = (solid lines) and of the same 
simulation with maximal convection (dashed lines). Top: entropy 
(red, left axis) and pressure (green, right axis). The entropy of the 
convective simulation has been raised by 0.5 dec to distinguish 
between the models. Bottom: density (red, left axis) and temper- 
ature (green, right axis). The positions of CFs in the adiabatic, 
non-convective simulations are marked with blue dashed vertical 
lines. 



the CF smears the discontinuity on a timescale that is set 
by the thermal velocity and the mean free path of the 
protons. The Spitzer m.f.p. A in an unmagnetized plasma 
with typical cluster densities is a few kpc, and depends on 
the thermal conditions on bot h sides of the discontinuity 
l|Markevitch fc Vikhlinir] 120071 ). Taking A ~ 10 kpc and a 
thermal velocity v ~ 1000 km sec -1 , and assuming that a 
CF is visible if it is sharper than L b s ~ 10 kpc, we get a 
characteristic timescale for CF dissipation, 



i -^obs olj obs , n 7 

t ~ — =— ~ ; — ~ 1U yr 



D 



(12) 



This result indicates that in order for shock induced CFs to 
be observable, either (i) they are formed frequently (e.g. , 
by a series of A GN bursts; see ICiotti fc Ostrikerl 120071 ; 
ICiotti et alj|2009h : or (ii) magnetic fields reduce the m.f.p. 
co nsiderably, as s ome evidence suggests (see the discussion 
in lLazarian l200fj ) . 

Heat flux driven buoyancy inst a bility (HBI; 
Parrish fc Quataertl l20Qg| ; iQuataerd 120081 : iParrish et al.l 



2009) tends to preferentially align magnetic fields perpen- 



dicular to the heat flow in regions where the temperature 
decreases in the direction of gravity. This effect has been 
argued to reduce radial diffusion within the inward cooling 
re gions in the cores o f cool core clusters. As suggested 
in IParrish fc Quataertl ()2008f ) . the temperature gradient 
across a CF is opposite to the direction of gravity, and HBI 
instability is expected to quickly form there (over 8 Myr 
in their example), aligning the magnetic fields parallel 
to the discontinuity and diminishing radial diffusion. 
This possibility needs to be addressed further by detailed 
numerical magneto-hydrodynamic simulations. Also, the 
timescale for achieving local thermal equilibrium between 
the electrons and ions below the virial shocks are long. It 
is unclear what the proper diffusion coefficients are in that 
case, and a particle transport analysis needs to be applied. 
These important issues are left for future work. 

CFs could also be degraded by subsequent shocks that 
sweep outwards across the CF. This is true regardless of 
their formation mechanisms. However, the CF contrast is 
only slightly diminished by this effect. This is shown ana- 
lytically by recalculating the Riemann problem described in 
§ [2] but starting from a general contact discontinuity with 
strength q instead of a primary shock. For q — 2 this decre- 
ment is between 4% for Mi < 2, to 20% for Mi > 10, 
as shown in fig. 1131 Such passing shocks seed Richtmyer- 
Meshkov instabilities which could cause the CF to break 
down. Although less efficient than Rayleigh- Taylor instabil- 
ities, these instabilities operate regardless of the alignment 
with gravity. The outcome depends on Mi and on any initial 
small perturbations in the CF surface. However, the insta- 
bility is suppressed if the CF becomes sufficiently smoothed. 

It is worth noting that KHI, which could break down 
CFs formed through ram pressure stripping and sloshing, 
does not play an important role in SICFs because there is 
little or no shear velocity across them. However, the sta- 
bilizing alignment of magnetic fields caused by su ch shear 
jMarkevitch fc Vikhlinir] 120071 ; iKeshet et ail I2009T ) is also 
not expected here. 
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Figure 13. Fractional decline (contours) in initial CF contrast q 
induced by a collision with a shock of Mach number M\ arriving 
from the dense CF side, for 7 = 5/3. The vertical dashed line 
corresponds to the maximal SICF contrast, q m ax- 



6.2 Observability 

lOwers et al] l|2009l ) present high quality Chandra observa- 
tions of 3 relaxed CFs. They all seem concentric with re- 
spect to the cluster centre, and spherical in appearance. 
They interpret these as evidence for sloshing, and find spi- 
ral characteristi c s in 2 of them. In the cold front sample of 
iGhizzardi et alj (|2010T ). 10 of the relaxed morphology clus- 
ters also host cold fronts. We argue here that some CFs of 
this type could originate from shock mergers. In the SICF 
model, CFs are spherical, unlike the truncated CFs formed 
by ram pressure stripping or sloshing, so no special projec- 
tion orientation is required. On the other hand, SICFs do 
not involve metallicity discontinuities, observed in some of 
these CFs. The o bserved contrasts of these three CFs in 
lOwers et al.l (|2009l ) are quite uniform, with best fit values in 
the range 2.0 — 2.1 (with ~ 20% uncertainty) in all three, as 
expected for SICFaJ. 

The following characteristics are specific to SICFs, and 
could thus be used to differentiate between SICFs and CFs 
formed by other mechanisms. These are also predictions for 
SICFs that are expected to form at radii that are not ob- 
servable today. 

Morphology. SICFs are quasi-spherical around the 
source of the shocks. In contrast, for example, a galactic 
merger defines an orbit plane, and the corresponding per- 
turbation (stripped material or sloshing and centre displace- 
ment) will cre ate CFs parallel to this plane. In some slosh- 
ing scenarios |Ascasibar &: Markevitcb] 120061 ) the CFs ex- 
tend considerably above the plane, but would never cover 

5 However, shocks sweeping over CFs could diminish their con- 
trast towards q ~ 2, regardless of CF origin. A larger, more com- 
plete sample of CFs is needed to determine the origin/s of the 
relaxed population. 



all viewing angles; an observed closed CF "ring" would 
strongly point towards an SICF. The statistical properties 
of a large CF sample could t hus be used to di s tingu ish be- 
tween the different scenarios. Ghizzardi et al. I |2010r i iden- 
tify 23 apparently relaxed clusters, 10 of which exhibit cold 
fronts. They find a perfect correlation between central en- 
tropy levels and cold fronts, with the 10 lowest central en- 
tropy clusters exhibiting cold fronts. This correlation is also 
con sistent with convection of low entropy gas from th e cen- 
tre (|Markevitch fc Vikhlini"nll2007l ; iKeshet et al.ll2009l ). but 
would require that all 10 CFs are viewed from a favorable 
viewing angle. We offer here an alternative interpretation, 
in which the low entropy and short cooling times imply rel- 
atively more active AGNs a nd consequently t he formation 
of shocks on shorter periods (|Ciotti et al.ll2009l ). The SICFs 
naturally cover 4-7T of the sky, alleviating the viewing angle 
statistical requirement. 

Amplitude. SICFs have distinct entropy and density 
contrasts that depend weakly on shock parameters; q is typ- 
ically larger than ~ 1.4 (assuming Mi ^ 2) and is always 
smaller than q ma x = 2.65. This is consistent with observa- 
tions. 

Extent. If cluster oscillations, as described in § 15.2 
are excited, SICF radii should be approximately logarith- 
mically spaced. Any shock that expands and collides with 
the virial shock will create an SICF at the location of the 
virial shock, far beyond the core. The external SICFs are 
younger, and would appear sharper owing to less degrada- 
tion due to the processes discussed in § 16.11 This concentric 
CF distribution thus resembles tree rings. Deep observations 
capable of detecting CFs at r ~ 1 Mpc are predicted to find 
SICFs (fig. I12J1 . Such distant CFs occur naturally only in 
the SICF model (observations by Suzaku may be able to 
probe this range wi th sufficient quality in the near future; 
iHoshino etalll2010h . 

Plasma diagnostics. Shocks are known to modify 
plasma properties in a non-linear manner, for example 
by accelerating particles to high energies and amplify- 
ing/generating magnetic fields. The plasmas on each side 
of an SICF may thus differ, being processed either by two 
shocks or by one, stronger shock. This may allow indirect 
detection of the CF, in particular if the two shocks were 
strong before the collision. For example, enhanced magnetic 
fields below the CF may be observable as excess synchrotron 
emission from radio relics that extend across the CF, in 
nearby clusters, using future high-resolution radio telescopes 
(MWA,LOFAR, SKA). 



7 SUMMARY AND DISCUSSION 

Our study consists of two parts. The first is a general, an- 
alytic discussion about cold fronts that form as a result of 
a merging between a primary and secondary shock prop- 
agating in the same direction. This is a novel mechanism 
to create cold fronts discussed here for the first time. The 
second part describes the possible relevance of this SICF 
mechanism in cluster environments using ID spherical hy- 
drodynamical simulations. 

We have shown in § [5] that when shocks moving in the 
same direction merge they generate a CF. The density con- 
trast across the CF is calculated as a function of the Mach 
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numbers of the two shocks. It is typically larger than 1.4 (if 
M > 2), and is always smaller than q ma x = 2.65. We sup- 
port the analytical calculation by a detailed investigation of 
a shock tube planar hydrodynamical simulation. 

In §[4] and §[5j using a ID spherical hydrodynamic code, 
we demonstrate that SICFs in clusters are a natural conse- 
quence of shocks that are generated at centres of clusters. 
We entertain two ways to invoke shocks: by injecting large 
amounts of energy near the centre (corresponding to AGN 
outbursts), and by abruptly changing potential the well of 
the cluster (corresponding to merger events). Finally, by 
simulating cluster evolution from initial cosmological per- 
turbations over a Hubble time we show that outgoing shocks 
that merge with the virial shock create very distinctive CFs. 
These shocks can be caused by a critical event (merger or 
explosion) but can also be invoked by stochastic oscillations 
of the cluster's core, caused by accretion of low angular mo- 
mentum substructure that perturbs the core. We show that a 
reverberation mode exists in the haloes of galaxies and clus- 
ters that causes periodic merging between the virial shocks 
and secondary shocks, producing an SICF every cycle. The 
simulated SICF contrast is consistent with the theoretical 
predictions of §[2] A more thorough investigation of this po- 
tentially important mode, including its stability in 3D is left 
for future works. 

We then discuss in § 16.11 the survivability and degra- 
dation in time of CFs after they are formed. The CF dis- 
continuity is smeared over time by diffusion, at a rate that 
depends on the unknown nature and amplitude of mag- 
netic fields. CFs are susceptible to heat -flux-driven buoyancy 
instability l|Parrish fc Quataertl 120081 . HBI; ), which could 
align the magnetic field tangent to the CF and potentially 
moderate further diffusion. SICFs, like all other CFs, are 
subject to Richtmyer-Meshkov instabilities from subsequent 
shocks passing through the cluster. Such collisions also re- 
duce the CF contrast until it reaches q ~ 2. Unlike most 
other CF models, an SICF is not expected to suffer from 
KHI. 

The predicted properties of SICFs are presented in 
i |6.2l and re p roduc e so me of the CF fea t ures discussed in 
lOwers et all |2009l ) and IChizzardi et all |2010l ). In partic- 
ular, we suggest that CFs in relaxed clusters, with no ev- 
idence of mergers, shear, or chemical discontinuities, may 
have formed by shock collisions. We list the properties of 
SICFs that could distinguish them from CFs formed by 
other mechanisms. The SICF model predicts quasi-spherical 
CFs which are concentric about the cluster centre, with con- 
trast q ~ 2, and possibly extending as far out as the virial 
shock. An observed closed (circular /oval) CF could only be 
an SICF. In the specific case of cluster reverberation, a dis- 
tinct spacing pattern between CFs is expected. It may be 
possible to detect them indirectly, for example as disconti- 
nuities superimposed on peripheral radio emission. 

Shocks originating from the cluster centre naturally oc- 
cur in feedback models that are invoked to solve the over- 
cooling problem. They are also formed by mergers of sub- 
structures with the BCG. Thus, SICFs should be a natural 
phenomenon in clusters. Further work is needed to assess 
how common SICFs are with respect to other types of CFs, 
and to characterize inner SICFs that could result, for exam- 
ple, from mergers between offset AGN shocks. The proper- 
ties of SICFs in 3D will be pursued in future work. 
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APPENDIX A: SICF PARAMETERS 

Using the thermal jump conditions, we find the density /temperature contrast across the CF given by 

P3l _ MlMl (7 + 1) [M 2 (7 - 1) + 2] / [(2M 2 - 1) 7 + 1] [(2Af 2 - 1) 7 + 1] ' 
9 " P3o M 2 [M 2 ( 7 - 1) + 2][Af?( 7 - 1) + 2] \ (7 + 1)[(2M 2 - 1) 7 + 1] 

where Mf ~ MqM\ is the Mach number of the newly formed shock. It is determined by including the velocity jump conditions; 
for 7 = 5/3 this yields 

1/5 = 1 Ml - 1 _ 4M 1 (M f - M )(M M f + 1) 

v/(5Af2 -l)(Af 2 + 3) Mf \J (5Mfi - 1)(M, 2 + 3)(5Af 2 - l)(Af 2 + 3) ' 



5M? 

4 L 



(5Af 2 - 1)(5M 2 - 1) 



which has three roots for Mf, only one of which is real. 

In the strong primary shock regime Mo S> 1, appropriate for example for a secondary merging with a virial shock, q is 
the root of 



q^+ql^ \^ = A + (Mo 2 ) , (A2) 
where 



.4 



7 + 1 \ ^ , ( 7+1 \^( Ml ~w) vwt + y~ 



7 - 1 + 2M- 2 J V7 - 1 + 2Aff 2 J V( 2M ? -1)7 + 1 



The contrast reaches its maximal value q ma x when M\ satisfies dA/dM\ — 0. This maximal value is a function of 7, as shown 
in Figure [3] 

For 7 = 5/3, Eqs. (| A2I simplifies to 



+-^-* 



4M 2 V V V 57 ^/(5M? -1)(M? + 3) 

Maximizing the contrast yields q max = 2.653, corresponding to Mi = 6.654. 

It is interesting to point out that the observed contrast q across an SICF could be used to impose an upper limit on 7, 
and hence constrain the equation of state. Figure [3] shows that q monotonically decreases with increasing 7 in the relevant, 
7 > 1.15 regime. After an SICF is produced with some initial contrast q ^ q m ax, its contrast typically declines in time as the 
CF is degraded by diffusion, convection, interactions with subsequent outgoing shocks, etc. Therefore, q m ax(si) > q constrains 
7. For example, an observed SICF contrast q > 2.653 (q > 3.108) would imply that 7 < 5/3 (7 < 4/3), suggesting the 
abundance of cosmic rays where the shocks collided. 



